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ABSTRA.CT 

A new  approach  to  the  formulation  and  analysis  of  stochastic 
models  of  chemical  reactions  is  presented.  Unim.clecular,  bimo- 
lecular,  and  enzyme  kinetic  reactions  are  considered  in  the 
irreversible  and  reversible  cases.  The  .methodology  is  based  on 
diffusio.n  apprc.x;i.T:ation3  and  represents  the  ti.me  evolution  of  the 
reaction  as  the  S'um  of  a deter.m.i.nistic  f'ur.ct  ion  and  a.n  Cr.nstein- 
Uhle.nbeck  process.  As  a result  the  marginal  dist ributions  are 
approxi.T.ately  Gaussia.n  'with  relatively  simple  mean  and  covariance 
para-meters,  and  the  dy.namic  behave  is  completely  characterized. 
The  stochastic  approach  which  uses  stochastic  differential  equa- 
tio.ns  is  a .'.atural  generalization  of  the  dete icninistic  approach 
which  uses  ordi.nary  differential  equations. 
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Ir.t  ro  Juct  ion 


In  this  paper  a new  approach  to  the  fomulation  and  analysis 
of  stochastic  models  of  chemical  reactions  is  presented.  Uni- 
molecular,  bimclecular,  and  enzyme  reactions  are  considered  in 
both  the  irreversible  and  reversible  cases.  The  reaction  types 
chose.n  Tor  a.nalysis  are  of  a si.mple  sort;  however,  this  is  for 
purposes  of  e;oposition  only,  and  the  methodology  can  be  easily 
exte.nded  to  very  co.mplicated  reactions.  The  results  developed 
i.n  this  paper  should  be  compared  for  their  simplicity  with  those 
obtai.ned  by  other  m.ethods.  The  survey  of  McQuarrie 
provides  an  excellent  summary  of  such  results. 

Early  development  of  the  stochastic  approach  to  chemical 
kinetics  is  due  to  Singer^^^  and  Bharucha-Reid^^^ . This  approach 
to  chemical  kinetics  typically  treats  the  concentratlo.ns  of 
various  molecule  types  as  a continuous  time  Markov  chain  with 
specified  tra.nsition  rates.  The  transitions  rates  are  used  to 
construct  the  Kolmogorov  forward  equations  which  characterize 
the  probabilistic  behavior  of  the  reaction.  The  exact  solution 
of  these  partial  differential  difference  equations  is  intractable 
except  in  the  simplest  case,  so  most  authors  are  content  to  find 
the  first  two  moments  of  the  process.  Even  this  often  leads  to 
differential  equations  whose  solutions  e.ntail  special  mathematical 
functions.  In  such  cases  it  is  difficult  to  gain  an  intuitive 
understanding  of  the  process.  Furthermore,  such  an  analysis  is 
very  limited  in  that  it  provides  information  only  about  moments 
and  not  the  exact  distribution.  Moreover,  it  provides  no  in- 
sight into  the  dynamic  behavior  or  n-dimensional  distributions 
over  time.  Finally,  as  the  reaction  becomes  more  complicated 
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the  analyols  becomes  more  intractable  and  the  Insights  less 
available . 

The  approach  presented  in  this  paper  is  approximate  rather 
than  exact,  but  overcomes  all  of  the  difficulties  which  arise 
with  a stochastic  analysis  mentioned  before.  Furthermore  it 
provides  a clear  connection  betv/een  the  usual  deterministic 
a.nalysis  and  a stochastic  analysis.  The  method  is  called  a 
diffusion  app roxi.mation,  and  examples  of  its  use  are  given  i.n 

ri^]  r^] 

the  worh  of  Gaver  ana  Lehoczky - ■'  and  McNeil  and 
Schach^^^.  The  approach  is  asymptotic  in  that  it  assiumes  the 
n.anber  of  molecules  involved  in  the  reaction  is  large,  although 
i.n  the  work  of  Gaver  and  Lehoczky  previously  cited  numerical 
studies  indicate  excellent  accuracy  in  small  systems  as  well. 

T.he  fundamental  idea  is  based  on  the  replacement  of  the  discrete 
state  space  by  a continuous  state  space  and  the  Markov  chain  by 
a Markov  process  (in  these  cases  a diffusion  process).  The  for- 
ward equatio.ns  (in  the  Markov  process  case  known  as  the  Fokker- 
Flanck  equations)  are  then  formulated  as  Ito  type  stochastic 
differential  equations  rather  than  their  equivalent  version  as 
partial  differential  equations.  This  formulation  permits  the 

use  of  the  Ito  calculus  for  the  manipulation  of  stochastic 

rgl 

differential  equtions.  The  reader  may  consult  Arnold'--^-'  or 
Gikhman  and  Skorokhod^^*^^  for  a treatment  of  this  theory.  The 
stochastic  process  is,  in  the  asymptotic  case,  transformed  into 
the  S'um  of  a deterministic  process  ( invariably  the  same  as  the 
process  arising  from  a deterministic  analysis  of  the  system) 
and  a stochastic  process  (Invariably  an  Ornstein-Uhlenbeck  process). 


% 
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This  characterization  of  the  system  provides  a full  description 
of  the  marginal  distribution  (Gaussian  with  specified  mean  and  c 
variance  functioiis)  over  time  and,  more  importantly,  of  the  n- 
dimensional  distributions  over  time  (Gaussian  with  specified 
mean  vector  and  covariance  structure).  The  characterization 
of  the  n-dimens ional  distribution  functions  is  important  for 
the  development  of  a proper  statistical  analysis  of  data 
gathered  from  such  a process  say  to  estimate  the  rate  constants 
of  the  reaction.  Since  the  stochastic  process  is  added  to  the 
deterministic  process,  the  deterministi c process  may  be  thought 
of  as  a first  order  approximation. 

It  is  important  to  note  that  the  methodology  given  in 
this  paper  can  easily  be  applied  to  very  complicated  reactions 
i.ncluding  those  with  many  compounds  and  steps.  The  results  are 
relatively  simple  and  easily  understood.  Once  the  nature  of  the 
Ornstein-Uhlenbeck  process  is  understood,  its  appeara.nce  in 
the  description  of  many  difference  reactions  allows  for  greater 
insig.ht  into  chemical  processes. 


2.  Unimolecular  Reactions 


In  this 
cons  lie  red. 

Me  Q earrie^^^ 
the  diffusion 
A.  Reaction 
ir  A(t) 
at  tirr.e  t, 
continuous  ti 


section,  two  si.Tiple  unimolecular  reactions  are 
The  two  selected  are  taken  from  the  review  by 
in  order  to  provide  a basis  of  comparison  for 
approximation  results  with  the  exact  results. 

represents  the  number  of  A molecules  presen 
then  [A(t),t<0]  is  typically  treated  as  a 
me  Markov  chain  with  transitions  given  by 


t 


Transition 


Rate 


a ^ a - 1 
a a 

all  others 


kadt  + o(dt) 

1 - kadt  +•  c(dt)  2.A.1 

o(  dt) 


where  o(h)/h  ^ o as  h ->  o. 

For  large  values  of  A(o),  the  time  between  tra.nsitions 
is  very  small  (having  a mean  value  of  k/A(t)).  Consequently, 
over  a short  period  of  time  many  transitions  will  occur.  It 
is  therefore  reasonable  to  assume  A(t  + dt)  - A(t)  = dA(t)  has 
appro.ximately  a Gaussta.n  distribution  with  mean  E(dA(t))  = 
-kadt  + o(dt)  and  variance  Var(dA(t))  = kadt  + o(dt).  One 
may  therefore  represent  A(t)  as  the  solution  of  the  Ito  sto- 
chastic differential  equation  (s.d.e.  in  the  sequel) 

dA(t)  = -kA(t)dt  + (kA(t))^/^  dW(t)  2. A. 2 


where  [W(t),  t > 0}  is  a standard  Wiener  process  and  dW(t)  = 
W(  t dt)  - W(t) . 
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The  state  space  for  the  Markov  chain  was  f 0, 1 , . . . , A( o) ] 
originally  and  is  now  [0,A(o)].  The  infinitesimal  mean  anh 
variance  in  (2. A. 2),  -kA(t)dt  and  kA(t)dt,  were  chosen  to 
match  that  of  the  original  Markov  chain  keeping  terms  of  order 
at  only.  The  approximation  of  the  discrete  process  by  a con- 
tinuous process  is  made  rigorous  in  the  work  of  Kurtz and 

ri6i 

Barbour'  . The  increments  are  modelled  as  having  a Gaussian 
aistributicn  in  view  of  the  central  limit  theorem.  I.n  the  spirit 
of  the  ce.ntral  limit  theore.m,  it  is  natural  to  introduce  another 


stochastic  pr'~ 
with  N = 

The  X.j  p 


[X.^(t),  t>0"(  where 
■ a(t)  is  an  arbitrary 
will  also  be  a diffusion 


s.a.e.  which  can  be  determined  using  Ito's 
[10,p.27]^  Specifically, 


X.,(t)  = 'A(t)  - :ia(t)/:i^/^ 
aetermi.nistic  f^unction. 


and  wil 


Lemma 


1 satisfy  a 
[9, P-70]  or 


dXjj(t)  = -N^^^(a(t)  + ka(t))dt  - kXj^.(t)dt 
-t-  (ka(t)  + 0(N'^'  “ ) ) 1W(  t) . 

As  A(o)  = N ^ 00  , the  term  in  (2.  A.  3)  will 
unless  its  coefficient  is  exactly  0.  The  function  a(t) 
therefore,  be  chosen  to  satisfy  the  equatio.n 


2.  A.y 


explode 


mu  3 1 , 


a'(t)  = -ka(t),  a(o)  = 1.  2. A. 4 

Given  that  a(t)  satisfies  (2.A.^)  and  therefore  equals 
a(t)  = exp(kt),  one  may  let  N -»■  oo  in  (2.A.]5)  and  conclude  from 
Barbour^^^^  Theorem  K that  [Xj^(t),  t>0]  converges  weakly  to 
{X(t),  t>0]  where  X satisfies  the  s.d.e. 


(2. A. 5) 


dX(  t) 
X(0) 


-k:X(t)dt  + 

0. 


The  lit’fuEion  approximation  consists  of  writing 

A(t)  a Na(t)  + ?I^'^^X.j(t)  :b  Ma(t)  + N^/^X(t).  {2.k.6) 

The  stochastic  process  X is  a nonstationary  Ornstein- 
Uhlenbeck  process  with  mean  C.  Many  facts  about  X can  be 
deduced  frcm  this  characterisation  including  the  Joint  distri- 
bution  functions,  the  moments,  and  the  spectr ‘ 

The  X process  will  have  mean  0 and  variance  given  to 

be  the  solution  of 

= -2kZ^  + kexp'-kt),  Z^,  = 0.  (2.A.?) 

Equation  (2. A. 7)  is  readily  solved  to  give 

Z^  = exp(-kt)(l  - exp(-kt)).  (2. A. 3) 

The  X process  is  Gaussian  and  X(t)  has  a Gaussian  distri- 
bution with  mean  expf-kt)  and  variance  exp( -kt) (1  - exp( -kt) ) , 
thus  A(t)  will  have  approximately  a Gaussian  distribution  with 
mean  A(o)exp(-kt)  and  variance  A( o) exp ( -kt) ( 1 - exp( -kt) ) . 

This  approximation  is  clearly  an  accordance  with  the  exact  results 
for  this  process,  namely  a Binomial  ( A( o) , exp( -kt) ) distribution, 
as  the  central  limit  theorem  guarantees  the  goodness  of  our 
approximation  even  for  small  values  of  A(o).  The  difi'usion 
approximat J on  (2. A. 6)  provides  a complete  description  of  the 
dynamic  behavior  of  the  reaction  rather  than  merely  describing 
the  static  or  marginal  behavior.  The  approximation  is  motivated 
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fcy  the  large  A(o)  case,  but  the  reader  is  referred  to  the  work 
of  Gaver  and  Lehocsky  ^ ^ ^ ^ ^ for  more  information  concerning 

the  accuracy  of  the  method.  Excellent  accuracy  has  been  reported 
for  values  of  A(o)  as  sm.all  as  20. 

?.  Reversible  Unimclecular  Reaction  A ^ ^ > 

[17] , [IS: , [19] _ ^2 

The  irreversible  reaction  of  section  2A  is  now  generalized 
to  the  reversible  case.  The  forward  and  backward  rate  constants 
are  k.,  and  ^2  respectively.  The  transitions  and  transition 
rates  appropriate  for  such  reactions  are  given  by 


Transit! on 
a ->  a + 1 

a -►  a - 1 

a -*•  a 

all  others 


Rate 

k2(B(o)  + A(o)  -a)dt+  o(dt) 

k^  adt  + o(dt) 

1 - (k^a+k^fBfo)  +A(o)  - a))dt  + o(dt) 
o ( dt ) 

2.3. 1 


The  A process  can  again  be  approxim.ated  by  the  s.d.e. 
aA(t)  = k2(B(o)  + A(o)  - A(t))dt  - k^A(t)dt  2.B.2 

+ (k2(B(o)  + A(o)  - A(t))^/^dW^(t)  - (k^A(t))^/^dW2{t) 


where  W^(t)  and  W2(t)  are  independent  standard  Wiener 
processes.  An  Xj^  process  is  again  introduced  where  = 

(A(t)  - Na(t))/N^'^^  and  N = A(o)  + B(o).  Let  = A(o)/N 

and  = B(o)/N.  Using  Ito' s Lemma  one  may  derive  the  s.d.e. 
for  the  process  to  be 


j 


2.  5.3 


3 


dX,j(t)  = (t)  - k2(l  - a(t)  ) + X^a(t)) 


- (k^  + k2)X.j(t)dt+  (^^(l  - a(t))  +k^a(t)  +0(rJ--/^))  iW't) 

1 /2 

Letting  N - « , the  coefficient  of  the  N*  term  must 
be  0,  so  a(t)  must  satisfy 


L/2 


,'(t)  = -k^a(t)  +k2(l  - a(t)),  a(o)  = a^ 


2.B.  4 


or 


4-  T- 


rC.  * -C, 


If  a(t)  is  given  as  above,  then  (X,j(t),  t>c]  converges 
weakly  to  {X(t),  t>0},  a diffusion  process  which  satisfies 


dX(t)  = -(k^  ^ k2)X(t)dt  + (k2(l  - a(t))+ k^a(t))^/^dW(t) 

x:o)  = 0. 


2.E.5 


The  X process  is  a nonstationary  Ornstein-Uhlenbeck 
process,  and  the  diffusion  approximation  methodology  treats  A(t) 
as  being  given  by  (A(o)  + B(o))a(t)  + (A(o)  + B(o))X(t).  One 
may  easily  compute  the  variance  of  the  approximation  process, 

Z^,  to  be  the  solution  of 

it  = -2(k^  + k2)Z^  + k2(l-a(t))  + k^alt"),  Z^  = C.  2.B.6 


The  solution  is  given  by 


Z^  = exp( -Kt) (1  - exp( -Kt) ) 


+ (k^  - k2)fa^  - k2/K) 


( exp ( Kt ) + 1 ) 


2.  B.7 


where  K = k^  + k2* 
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I One  may  let  t -»■  « to  deluce  the  equilibrium  ■list ribution 

of  A(t).  In  thio  case  A(t)  will  have  a steady  state  Gaussian 
distribution  'with  mean  k2(A(o)  + 3(o))/K  and  variance 
k.ik.2(A(o)  + B(o))/k'^.  These  results  extend  those  previously 
obtained  for  this  reversible  reaction.  Of  particular  importa.nce 
is  the  characterization  of  the  dynamic  behavior  in  terms  of  an 
Ornste in-Uhle.nbeck  process.  The  n-dimensional  distribution 
f.uictions  over  time  are  thus  given  by  an  n-dimensional  Gaussisun 
distribution. 

Many  other  unimclecular  reactions  can  be  treated  in  this 
man.ner  including  the  parallel  first  order  and  triangular  reac- 

m 

tions  described  by  McQuarrie'-  . Instead  we  next  treat  the  .more 
complicated  nonlinear  situations. 
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^ . Bimolecular  Reactions 

The  bimolecular  reactions  to  be  considered  are  again  siummarized 

n 1 

in  the  survey  of  McQuarrie*-  . The  methodology  of  diffusion  approx.i- 
m.ations  is  especially  useful,  because  it  provides  sLmple  characteri- 
zations of  the  evolution  of  these  processes  in  terms  of  the  sum  of 
a deterministic  function  and  an  Ornste in-Uhlenbech  process.  This 
is  in  sharp  contrast  to  the  results  obtained  from  a.n  exact  analysis. 
Such  results,  when  they  can  be  obtained  at  all,  i.nvariably  involve 
complicated  ccmbinations  of  special  mathematical  functions.  Con- 
seque.ntly,  it  .has  bee.n  i.mpossible  to  obtain  any  i.ntuitive  u.nder- 
standing  of  the  stochastic  nature  of  the  che.mical  reactio.n. 

’x 

A.  The  Irreversible  Reaction  2A  ->  B . 

A number  of  asymptotically  equivalent  stcc.hastic  models 
are  possible.  If  A(t)  represe.nts  the  concentration  of  A 
.molecules  at  ti.me  t,  then  [A{t),  t>0]  is  taken  to  be  a 
Markov  chain  with  tra.nsitions  a.nd  transition  rates  give.n  by 


Trsnsi  ‘ticn 


a a - 2 


a a 


all  others 


Rate 


ka"^ 


2(A(o)  + E(o)) 


dt  4-  o(dt) 


1 - 


, 2 
xa 


2(A(o)  + B(o)) 


dt  + o(dt) 


3.  A.l 


o(  dt) 


Using  (3. A.l)  and  defining  dA(t)  = A(t-t-dt)  - A(t),  it  is 
easy  to  show  E(dA(t))  = -kA^( t ) dt/( A( o ) + B(o))  + o(dt)  and 
2 

Var(dA(t)  = dt  + o(dt). 

A(o)  + B(o) 


U4 
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One  xay  therefore  approximate  the  behavior  of  the  A process  using 
the  s.a.e. 


,(t)  = - Miiil  It  + ..llj  iw(t) 


3.  A.  2 


N 


N 


v;here.  A(o)  + 5(o)  = N. 


Following  the  methods  of  section  2,  the  process  X.j(t)  = 
1 /2 

(A(t)  - Na(t))/N  is  introduced,  and  using  Ito' s Lemma  the 

s.d.e.  governing  [X^,(t),  t>0]  is  found  to  be 

dX,,(t)  = -N^'^^(a'(t)  + ka^(t))dt  - 2k.a(  t)X^j(  t)dt 


.2  / 4.  s , 1/2 


+ f2ka^(t)  + 0(N‘ 


1/2 

))  dW(t). 


X A 


' y2 

The  coefficient  of  must  be  identically  0,  thus  a(t) 

must  satisfy 


^.A.l 


a(  o) 


a'  (t)  = -ka  (t) , a(o)  = A(o)/(  A(c)  + 3(o) ) or  a(t)  = i + 

Assuming  a(t)  satisfies  (3-A.'^)  and  N ^ «>,  [ X^,j(  t) , t > 0} 

converges  weakly  to  [X(t),t£_0)  which  satisfies 


dX(t)  = -2ka(t)X(t)  + (2ka^(  t)  )^/’^dW(  t) 
X(o)  = 0. 


^ • A.  5 


Once  again  the  diffusion  approximation  gives  a representation 
for  A(t)  as  (A(o)  + B(o) ) ( a( t)  + X( t) ) where  a(t)  is  given 
by  (j.A.'^i).  The  X process  Is  a nonstationary  Ornstein-Uhlenbeck 
process  with  variance  satisfying 

Z^  = -Uka(t)Z^  + 2ka^(t),  Z^  = 


0. 


J.  A.  6 
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I Equation  (3- A. 6)  is  easily  solved  to  give 

Z = f a(o) ( (1 + a(o)kt)^  -1 )/(l  + a kt)^.  3-A.7 

It  follows  that  marginally  A(t)  v/ill  have  approximately  a 
normal  distribution  with  mean  (A(o)  + E(o))a(t)  and  variance 
'A(o)  •<'  B(o))Z*  whei'e  a(t)  and  Z.  are  given  in  (3.A.Z)  and 
(3. A. 7). 

The  dynamic  behavior  of  the  A process  can  be  easily  deduced 
from  its  characterization  as  an  Ornste in-Uhlenbech  process.  The 
sim.plicity  of  this  characterization  should  be  co.ntrasted  with  the 

n 1 

analysis  presented  by  .McQuarrie'-  . Only  the  moments  are  presented, 
and  they  appear  as  co.mplicated  cc.mbinations  of  gamma  functions. 

It  should  be  .noted  that  the  stochastic  formulatio.n  give.n  i.n 
(3.A.I)  differs  fro.m  that  give.n  by  all  other  researchers  with  the 
prese.nce  of  the  A(o)+B(o)  factor  i.n  the  denominator  of  the 
tra.nsition  rates.  This  factor  is  important  when  deali.ng  with  non- 
linear reactions.  First,  this  denominator  .makes  the  units  of  the 
I rate  reaction  consta.nts  the  same  in  the  unimolecular  a.nd  bim.olecular 

cases,  whereas  they  would  be  different  if  the  transition  rates 
were  taken  to  be  kA  (t)dt.  Second,  it  must  be  the  case  that  as 
the  reaction  proceeds,  B molecules  are  produced  and  this  makes  it 
difficult  for  the  A molecules  to  be  paired  up.  Using  the  model 
without  the  A(c)  + B(o)  factor,  the  time  for  this  reaction  to  be 
50'  complete,  T^,  is  l/2A(o).  If  the  initial  number  of  A 
I molecules  is  doubled,  the  time  for  this  reaction  to  be  75'  co.m- 

t 

I Plete,  3/^A(o).  In  fact  the  latter  time  should  be  larger 

than  the  former,  since  for  the  latter  to  be  75'?  complete,  it  must 

I i 
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first  be  50-  complete  (at  which  time  there  are  A(o)  A molecules 
left).  This  tak.es  time  T^  = 1/m-A(o).  After  reachi.cg  50^  com- 
pletion it  takes  an  additional  T,^  time  lunits  to  reach  75''- 
ccmDletio.c.  Unfortu.nately  T, , = T.  so  should  exceed  T., 

which  leads  to  a clear  contradiction.  Throughout  this  paper,  we 
.normalize  all  nonlinear  transition  rates  by  the  number  of  .molecules 
present.  This  will  equalize  the  transition  rates  whe.n  both  linear 
and  nonlinear  terms  are  included  (see  section  30). 

5.  The  Irreversible  Reaction  A+B^  g[l]  > [13] , ^20] 

The  quantity  A(t)  again  represents  the  number  of  A 
molecules  at  time  t.  We  let  A(o)  + B(o)  = N a.nd  note  that 
B(t)  = B(o)  - A(o)  + A;t)  represents  the  number  of  E .molecules 
prese.nt  at  time  t.  The  transitions  for  the  Markov  chain  (A(t), 
t>0]  are  given  by 


Transitions 
a -*  a - 1 
a -*■  a 

all  others 


Rate 

ka(B(o)  - A(c)  -h  a)dt/:i  + o(dt) 

1 - ka(B(o)  - A(o)  -f  a)dt/N  -f  o(dt) 
o(dt) . 

3.  B.l 


The  approximate  s.d.e.  for  the  A process  is  given  by 

1 /2 

dA(t)  = (B(o)  - A(o)  + A(t))dt  + (k^^ip-^  (B(o)  - A(o)  + A(t)))  dW(t) 

3.B.2 

1/2 

The  process  Xj^(t)  = (A(t)  -Na(t))/N  ' will  also  be  a 
diffusion  process,  and  using  Ito' s Lemma  the  s.d.e.  for  the 
[Xj|(t),  t^O")  process  is  found  to  be 


Ih 

+ ka(t)(2b  - 1 - a(t)))dt 

- k(2b  - 1 + 2a(  t)  )X^j(t)dt  J.E.J 

1/2 

+ ( ka(  t ) ( 2b  - 1 + a(  t ) ) + 0(X”^'^^))  Ik/t) 


where  b = 5(o)/(A(o)  + 5(o)). 

Letting  X <»  and  assuming  a(t)  satieties 

a'(t)  = ->a(t)(2b-l  + a(t)),  a(o)=l-b  J.3.1 

then  (X,j(t),  t>0  ) converges  weakly  to  {X(t),  t>0}  v/hich 
satieties 

- 

iX(t)  = -k( 2b  - 1 + 2a( t) ) X( t) dt  + ( ka( t ) ( 2b  - 1 + a( t ) ) )‘  "dw(t) 
X(0)  = C.  1.B.5 


The  deterministic  function 

a(  t) 

is  easily  fO-;nd 

to  be 

(2b  - 1) (1 



— , 0 < b <;  1,  b 

/ 1/2 

a(  t)  = < 

bexp ( (2b  - 1 ) kt) 

« ' -1 

b) 

3.B 

1/(2  + kt) 

, b = 1/2 

The  X process  is  nonstationary  Crnstein-Uhlenbeck,  hence 
the  dynamic  behavior  can  be  easily  deducea.  The  mean  will  be 
0 tor  all  t,  and  the  variance,  will  satisfy 

Kj 

= -2k(2b-l  + 2a(t))Z^  + ka(t)(2b-l  -f  a(t))  3-E.7 

Z^  = 0 . 

This  equation  can  be  easily  solved  using  (3-B.o)  to  give 
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a^t) 


.3.3 


[a(t)Oxa(t))l 


,2 


0 -r  2 a (’  t ^ 


0 + 2a(o) 


L<3‘^a(  t ) ( 0 + a(  t ) ) 


- — rlog 


/ S + a( t) 

a(ol 

1 a(t) 

0 - a(o)  j 

0‘^a(o)  ( 0 -r  a''o) ) 

0 < b < 1 


0 = 2b  - 1 

The  stcchastia  evoluticr  oT  A is  thus  give.a  in  terms  of 
a.n  Orr.ste in-Uhls.nbe  ck  process.  For  large  values  of  A(o)  and 
3(o),  A.t)  will  be  appro .xiinately  normally  distributed  witn 
mea.n  (A(o)  + 3(o))a(t)  and  variance  (A(o)  + E(o))Z^  where 
a(t)  a.nd  are  give.n  above.  The  si.mplicity  of  these  results 

sho.;ld  be  co.ntrasted  with  the  complexity  of  the  exact  analysis. 

C.  Reversible  Bio.molecular  Reactions^^^ , [15] , [IS] , [21] , [22] , [23] , [24] 
The  results  given  in  sections  3A  and  33  can  be  easily 
extended  to  tne  case  where  the  reaction  is  reversible.  In  the 
case  of  a reversible  reactio.n  one  may  also  calculate  an  equili- 
brium or  steady  state  distribution.  As  usual  A(t)  represents 
the  concentration  of  A molecules  present  at  ti.me  t,  and 
the  concentration  of  all  other  types  of  molecules  can  be  deter- 
mined from  A(t)  and  the  initial  conditions.  We  consider  four 
types  of  reversible  reactions 


I.  A + R « c D 

’-^2 


II.  A + 3 . -4  C 


III.  2A 


2A, 


+ D 


^2 

k. 


IV. 


ic. 


The  method  ot’  diffusion  approximations  outlined  in  the 
previous  sections  can  be  routinely  applied  to  each  of  the  four 


reaction  types  given  above.  The  details  of  the  derivation  are 
omitted  as  they  are  identical  to  those  give.n  in  section  5A  snid 
3B.  It  is,  however,  i.mportant  to  .note  that  all  ncnli.near  tran- 
sition rates  will  be  normalized  by  N,  the  total  .number  of 
molecules  present  at  time  0.  This  is  especially  i.mportant  i.n 
cases  II  and  IV.  If,  in  II,  we  assume  A(o)  and  3'o) 
are  of  a comparable  order  of  magnitude,  the.n  the  u.n.nor.mal iced 
forward  tra.nsition  rate  will  be  an  order  of  magnitude  larger 
than  the  backward  rate.  Consequently,  the  reaction  will,  for 
large  values  of  A(c),  be  esse.ntially  equivalent  to  the  irre- 
versible reaction.  The  equilization  of  transition  rates  is 
acco.mplished  by  normalization. 

•‘■he  process  [A(t),  t>0]  is  taken  to  be  a Markov  process 

with  appropriate  transition  rates.  The  diffusion  approximation 

1 /2 

methodology’  treats  A(t)  = Na(t)  + N X(t)  where  N is  the 
number  of  molecules  of  all  ty'pes  at  time  0,  a(t)  is  a deter- 
ml.nistlc  function,  and  X(t)  is  a diffusion  process  satisfying 
a s.d.e.  We  ’«frite 


a'(t)  = 

f^a^(t)  - f 

2a(t)  + f^. 

a(o)  = a^ 

5. 

C.l 

dX(t)  = 

-g( t)X( t)dt 

+ (h(t))^/2 

dW(t) 

X(0)  = 

0, 

3- 

C.  2 

let  a^ 

II 

o 

O' 

o 

= B(o)/N,  c 

o = C(o)/N, 

and  d = D( 
o 

o)/N 

We  stummarize  below 

the  deflnisions  of  f^. 

^2’ 

and 

for  each  of  the  four  types  of  reactions. 


The  mean  of  A(t)  is  obtainable  as  :'Ia(t)  where  a(t) 
is  the  solution  of  (3.C.I).  The  steady  state  is  obtained  by 
solving  v/ith  a' ( t)  = 0 and  a(t)  = a.  The  variance 

of  Ait),  is  obtained  by  solving 

= -2g(t)Z^^  + h(t),  = 0.  5.C.3 

In  all  cases  A(t)  will  be  a Gaussian  process.  This 


i 


IS 


characterization  is  especially  appealing  when  compared  with  the 
exact  answers  where  the  moments  alone  require  complicated  combi- 
nations ol'  hypergeometric  functions. 

a.  Enzyme  Ki.net ic  Models 

The  last  class  of  ki.netic  models  to  be  addressed  are  those 
relating  to  enzyme  kinetics.  Two  simple  enzym.e  reactions  are 
studied,  one  which  leads  to  the  classical  Michaelis-Menten  or 
3riggs-Haldane  theory  and  a second  which  includes  enzyme  inhi- 
bitors. It  is  the  purpose  of  this  section  to  provide  a new 
stochastic  analysis  of  these  reactions  one  which  may  be  compared 
to  that  given  by  Bartholomay ^ and  Jachimowski  et 
The  stochastic  process  (now  multivariate)  will  again  be  a.n 
Ornstei.n-LTlilenbeck  process.  The  moments  of  the  reaction  and  the 
n-di.mensional  distributio.ns  will  ce  characterized  exactly. 


A.  3 + E; > ' SE] »P  -1-  E 

The  simple  e.nzyme  reaction  given  above  gives  rise  to  the 
Mic.haelis-Menten  or  Briggs-Halda.ne  theory  which  postulate  on 
intermediate  [SE]  complex.  Let  S(t)  represent  the  substrate 
concentration  at  time  t,  E(t)  represent  the  enzyme  concentra- 
tion at  time  t,  and  C(t)  and  P(t)  represent  the  concentra- 
tions of  enzyme-substrate  complex  and  products  at  time  t. 
{(S(t),  E(t)),  t>0]  is  taken  to  form  a Markov  chain  since 
C(t)  = E(o)  - E(t)  and  P( t)  = S(o)  - S{ t)  - C(  t) . 

Transition  Rate 

k^SE 


(S,E)  ^ (S-1,E-1) 


E(o) 


dt  + c(dt) 


lo 


Transition 

( continued) 

Rate 

( S-k1,E-1) 

k2(E(o) 

- E) it+  o( dt) 

(S,E+1) 

kj(E(o) 

- E) dt+  o( dt) 

( 3,  E) 

\E(o) 

(k2  + k,)(E(o) 

all  others 

o( dt) . 

4.  A.l 


E)  It  + c(lt) 


The  transitions  and  rates  given  in  f‘i.A.1)  cn  be  used 
to  give  the  s.d.e.  representation  of  the  reaction.  Specifically 


dS(t)  = -k  dt  + k^(E(o)  - E(t)  )dt  , , P 

^ E(c)  ^ ^ 

/k^S(t)E(t) .1/2 

dW^(t)  - k2(E(o)  - E(  t)) 

\ 2(o)  / / 


lE(t) 


k. S(t)E(t) 


+ k2(E(o)  - E(t))dt  + k^(E(c)  -E(t))dt 


/k,S(t)E(t)  I 

clW^(t)  + k2(Eo)-E(t))  ^/^dW2(t) 

\ E(o)  ' 

+ jk^(E(o)  - E(t))  j^/“dW-(t). 

We  next  introduce  Xj^(t)  = (S(t)  - Ns ( t)  and  Y,j(t)  = 

(E(t)  -Ne(t))/N^^^  where  N = E(o).  The  s.d.e.'s  governing, 
i(Xjj(t),  Yj,j(t)),  t>0]  can  be  easily  deduced  usi.ng  the  Le.tina  to  be 

dXj^ft)  = (t)  +k^s(t)e(t)  - k2;i  - e(t)))dt 


- k^(s(t)Y.^(t)  +e(t)Xj^(t))dt  - k2Y^(t)dt 

-1/2 

-(k^s(t)e(t)  +0(N  dW^(t) 

1/2 

+ (k2  (1  - e(t))  +0(N'^/^))  dW2(t) 


4.  A.  3 
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iY^(t)  = -N^^^(e'  (t)  +k^3(t)e(t)  - ( ^2  + k,  ) ( 1 - ef  t) ) ) dt 

- k^(s(t)Yj^(t)  + e(  t)X,j(t) ) It  - (k^  + k^)Y,j(t)dt 

•i  1/2 

- (k^a(t)e(t)  + 0(:r*''^^) ) dW,  (t) 

1/2 

+ (k2(l-e(t))  + 0(:r^/^))‘  dW2(t) 

1/2 

+ (k,(l-e(t))  + 0(rr^^^))  dW^(t). 

The  analysis  proceeds  by  letting  N -»  « and  requiring 
s(t)  and  e(t)  to  satisfy 

3'(t)  = -k.s(t)e(t)  + k2(l-e(t)) 

h.A.ii 

e'(t)  = -k^s(t)e(t)  ^ (k2  + k,)(l  - e(t)) 

efo)  = 1,  s(o)  = S(o)/E(o). 

Of  course  can  be  augmented  by  the  equation 

c'(t)  = k^e(t)s(t)  - (k2 + k^) (1  - e( t) ),  and  this  system  of  dif- 
ferential equations  is  now  a normalize i version  of  the  Michaelis- 
Menten  equtions.  The  Michaelis-Menten  theory  can  thus  be  thought 
of  as  a first  order  approximation  of  the  enzyme  reaction-  The 
second  order  approximatio.n  arises  from  an  analysis  of  the  sto- 
chastic term. 

If  s(t)  and  e(t)  satisfy  (k.A.U-),  then  as  N -* 

[(Xj^(t),  Yj^(t)),  t>0]  converges  weakly  to  [(X(t),  Y[t)),  t>0] 


which  is  a diffusion  given  by 


21 


A(t)  = 


k^e(t) 

kj^e(t) 


kj^s(t)  + 

(t)  + k2  + k^ 


B(t)  = 


(k^ s(t)e(t))^  (k2(l-e(t)) 

/'I. 


-(k^s(t)e(t))^  (k2(l-e(t)))^  (k3(l-e(t))) 


The  equation  given  by  (4. A. 5)  is  a nonstationary  bivariate  linear 
equation.  It  is  easy  to  show  that  the  eigenvalues  of  A(t)  have 
strictly  negative  real  parts,  hence  the  process  is  Ornstein-Uhlenbeck. 
The  marginal  distribution  of  (X(t),Y(t))  will  be  Gaussian  with  mean 
(0,0)  and  covariance  matrix  given  to  be  the  unique  nonnegative 

definite  solution  of  the  matrix  Riccatl  equation 


It  - * ?c5?-  lo  - S 


4.  A.  6 


The  above  equation  can  be  rewritten  as  follows.  Let  a^^Ct)  = 
Var(X(t)),  a2<t)  = Var(Y(t)),  and  0^2^^^)  “ Cov(X(t)  , Y(t)).  Then 
calling  V(t)  = (0^ (t)  , a^2 find 


V(t)  = -pj(t)  + 


V(0)  = 0 


4. A. 7 


where 


and 


?t 


2kj^e(t) 

k^e(t) 


2(k^s(t)  + k2) 


kj^(e(t)  + s(t))  + k2  + k^  k^s(t)  + k2 


2k^e(t) 


'k^s(t)e(t)  + k2(l-e(t)) 
k^s(t)e(t)  + k2(l-e(t)) 
kj^s(t)e(t)  + (k2  + k^Kl-eCt)) 


/ 


2(kj^s(t)  + k2  + k3> 


/ 
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Equation  (4. A. 7)  may  be  solved  formally  to  give 
V(t)  = exp(-  p^du)  exp(F^/^  p^ds)du 


4. A. 8 


which  can  be  computed  using  standard  numerical  integration  techniques. 

The  diffusion  approximation  methodology  thus  provides  a complete 
description  of  the  enzyme  reaction  in  terms  of  a bivariate  Ornstein- 
Uhlenbeck  process.  It  appears  that  no  closed  form  expressions  are 
available,  so  s(t),  e(t),  and  V(t)  must  be  evaluated  using  numerical 
integration.  It  is  possible  to  obtain  approximations  by  adopting  the 
Briggs  and  Haldane  steady  state  assumption,  c'(t)  = 0,  This 
sinplif ication  allows  for  the  derivation  of  certain  closed  form 
expressions , 

It  should  be  noted  that  the  stochastic  description  of  this 
enzyme  reaction  is  valuable  not  only  for  the  insight  it  provides 
into  the  process,  but  also  it  is  necessary  for  the  development  of  a 
proper  statistical  analysis  of  reaction  data,  say  to  estimate  the 
Michaelis-Menten  constant  (k2+k.2)/k^.  The  analysis  presented  can  be 
easily  generalized  to  far  more  complicated  reactions  including 
reactions  with  several  intermediate  enzyme- substrate  complexes, 
several  substrates,  and  reversible  reactions. 


k^  k,  k, 

B.  Enzyme- Inhibitor  Reaction  [27]  E+S  — >[ES]— -»E+P  E+I  [ El  1 

^2  ^5 


We  conclude  with  a simple  enzyme-inhibitor  reaction.  The 
analysis  of  this  system  parallels  that  given  in  section  4A,  except 
that  a higher  dimensional  process  must  be  used.  Again  the  diffusion 


approximation  analysis  leads  to  a characterization  in  terms  of  a 
multivariate  Omstein-Uhlenbeck  process. 

Define  S(t),E(t),  and  I(t)  to  be  the  concentrations  of 
substrate,  enzyme,  and  inhibitor  at  time  t.  If  Cj^(t)  and  C2(t) 
represent  the  concentrations  of  [ES]  and  [El],  then  both  can  be 
calculated  from  the  initial  conditions  by  C2(t)  = I(o)  - I(t)  an 
C^(t)  = E(o)  - E(t)  - C2(t).  {(S(t),E(t),I(t)),  t>0}  can  be 

treated  as  a Markov  Chain  with  transitions  and  rates  given  by 
Transition  Rate 


(S,E,I) »(S-1,E-1,I) 


(S+1,E+1,I) 

(S.E+l.I) 


(3,E-1,I-1) 


(S,E+1, I+l) 


k.,  S Edt 

i-  o(dt) 

E(0) 

k,,(E(0)-E-C.,)dt  + c(dt) 

^ ^ 4.B.1 

k3(E(0)-E-C2)dt  + o(dt) 

k^EI 

— - — dt  + o(dt) 

E(0) 

k5(I(0)-I)dt  + o(dt) 


with  C2  = 1(0)  - I. 

We  omit  writing  the  approximate  s.d.e.'s  for 
(S(t),E(t),I(t)).  Let  Xj^(t)  = (5(t)  - Ns(t))/N^,  = (E(t)  - 

Ne(t))/N^,  and  Zj^(t)  = (I(t)  - Ni(t))/N^  where  N = E(0).  An 
asymptotic  analysis  indicates  that  the  deterministic  functions 
s(t),  e(t),  and  i(t)  must  satisfy 


s' (t)  - -k^  s(t)e(t)  + k2(i-e(t)  - iQ  + i(t)) 

e'  (t)  * -k^  s(t)  e(t)  + (k2  + k3)(i-e(t)  - i^  + i(t))  4.B.2 

i'(t)  - -kiie(t)i(t)  + k5  (i^  - i(t)) 


where  e(0)  = 1,  I(0)/E(0)  = , and  = S(0)/E(0). 

The  stochastic  noise  process  { (X^ ( t)  , ( t)  , ( t)  ) , t^O} 

will  converge  weakly  to  {U(t)  , t^O^  where  U(t)  = (X(t)  , Y(t)  ,Z(t) 
and  will  satisfy 

d U(t)  = -A^U(t)dt  + B^d  W(t),  U(o)  = 0 4.B.3 


with 

rt 

II 

(W^(t),W2(t).W3 

l(t),W4 

(t)  ,W3(i 

1 

' k^e(t) 

k^s(t)  + ^2 

k 

2 \ 

t > 
rr 

II 

k^e(t) 

k^s (t)  + k2 

+ k^  + 

k^i(t) 

-0<2 

+ k^  - k^) 

1 

1 0 

C)  + kj  j 

g2<t) 

0 

0 

?t  = 

g2(t) 

g3(t)  - 

g4(t) 

g5(t) 

0 

0 

0 

• g4<t) 

g5(t)  J 

with  gj^(t)  = (k^s(t)e(t))^,  = (k2(l-e(t)  - i^  + i(t)))^, 

g^Ct)  = (k^Cl-eCt)  - i^  + i(t)))^,  g^(t)  = (k^e(t)i(t))^,  and 
g5(t)  = (k3e(t)i(t))^. 

The  eigenvalues  of  have  strictly  negative  real  parts, 

hence  the  U(t)  process  is  0ms tein-Uhlenbeck  with  mean  0.  It 
follows  that  U(t)  will  have  a Gaussian  distribution.  The  covariance 
matrix  at  time  t,  is  given  as  the  unique  nonnegative  definite 

solution  of  (4. A. 6).  The  diffusion  approximation  then  gives  (S(t), 
E(t),  I(t))  to  be  aproximately  E(o) (s (t) , e (t) , i (t) ) + 

(E(o))^  (X(t)  , Y(t)  ,Z(t))  where  (s (t)  , e(t)  , i (t) ) is  given  by  (4.B.2), 
and  the  stochastic  part  is  specified  by  (4.B.3).  It  appears  impossible 
that  s (t) ,e (t) , i(t) , and  can  be  found  in  closed  form;  however, 

one  can  easily  use  standard  numerical  procedures  to  compute  them.  It 


25 


is  possible  that  further  approximations  similar  to  the  Briggs- 
Haldane  steady  state  assumption  can  be  invoked  to  obtain  some 
closed  form  information.  The  joint  distribution  of  the  process 
over  time  will  be  given  by  a trivariate  Gaussian  distribution. 
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20.  Abstract 


differential  equations  is  a natural  genera  I i za t i on  of  the  deterministic 
approach  which  uses  ordinary  differential  equations. 
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